Could you help me understand why the residue number should be hardcoded instead of looked up from the topology? Sometimes residue information is meaningless but it's often useful, and modifying it at export time is something I'd like to avoid

We cant' always start residue ID at 1, this is often a different value in user inputs and ought to be respected, i.e. this protein which starts a 3 for whatever reason https://github.com/openforcefield/protein-ligand-benchmark/blob/a0b67b93c4f12f63a8bb417a1f6e4df4dc10c6c5/data/shp2/01_protein/crd/protein.pdb#L591-L621

Believed you had moved it to #1016 or #1024 / planned on encountering it there. I'll review this again then.

It doesn't matter, but I would want to know "whatever reason" this example starts at 3 if this were the fix.

I don't know why it's the case, but it's common enough in biophysics that we need to support it, i.e.

• https://github.com/openforcefield/protein-ligand-benchmark/blob/a0b67b93c4f12f63a8bb417a1f6e4df4dc10c6c5/data/tyk2/01_protein/crd/protein.pdb#L332-L337
• https://github.com/openforcefield/protein-ligand-benchmark/blob/a0b67b93c4f12f63a8bb417a1f6e4df4dc10c6c5/data/cdk2/01_protein/crd/protein.pdb#L623-L628

PLB is a relatively high-quality dataset that's valued, if not endorsed in some capacity, by consortium alumni and industry partners. It confuses me, like plenty of PDB files, but it's at least a concrete starting point to work on.

Note that this applies to residue IDs defined in biopolymers, not necessarily applicable to non-biological polymers and certainly not applicable to cases in which Interchange applies "residue IDs" on the fly just to make GROMACS happy. One of the details here that's been a thorn in my side for a decade or so is that these engines inconsistently handle residue and molecule indices, almost none handling them both at the same time.

@mattwthompson
An additional loop appears to be necessary within the existing one which covers unique_molecules in order to loop through the topology_mol_idx in the lists:
example: 0: [. (0, {0: 0, 1: 1, 2: 2}), <--the existing loop only accesses the unique_molecule but we need the list int.
(1, {0: 1, 1: 0, 2: 2})]...
This can be achieved via:

        list_of_mol_tuples = unique_molecule_map[unique_molecule_index]
        for (topology_mol_idx, atom_map) in list_of_mol_tuples:
            unique_molecule_index = topology_mol_idx    # lazy diction assignment

...then in GROMACSAtom:
residue_index=unique_molecule_index + 1

This change prints the maximum molecule int in the system (originally all the minimum) which could indicate a few things:

1. Additional loop logic is needed.
2. The "writer" either only checks the final unique_molecule or is overwritten by each new identical molecule.
3. GROMACSMolecule contains one GROMACSAtom set/molecule instance that overwrites itself before sending to the system.

I checked the last possibility by reviewing the contents of GROMACSMolecule after each append below GROMACSAtom (I checked this a number of ways, but this seems to be the most telling. Printing the GROMACSMolecule ID is a little confusing):

            print(f"After append: unique_molecule.name = {unique_molecule.name}")
            for a in molecule.atoms:
                print(a)

This result, interestingly, does increment the residue_index correctly (testing with the original loop will only throw one instance of course) which suggests it may be the writer or a process nearby upstream. Regardless, I think the loop here does need to be changed, possibly in addition to a writer fix.

Also, both loops in _gromacs.py may need to be changed (line 113), but the errors I got from not doing so were inconsistent.

I think I'm quite close and I could try to track down the writer in comp/interchange which my debugger accesses as the next path, but you may be able to solve this far faster than me. Let me know your thoughts. I spent a while on this and to little avail since I accidently tracked deepcopy and add_mol_keep_cache for molecules before returning to this file thinking it would help me understand the background.

To be frank, there's too much going on here for me to follow, and I don't fully understand the behavior you're trying to fix. If you can open an issue with a minimally reproducing example demonstrating a difference between expectation and behavior, that would greatly help me understand the objective here.

I think #1024 there are changes you might be after? In particular I think this test looks at residue IDs in .gro files when multiple copies of a unique molecule are present

https://github.com/openforcefield/openff-interchange/pull/1024/files#diff-796291360ccaaab03d27aab793674a9c3b447a64c3f6fd2402b600611f1e36f5R200-R214